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ome  rules  for  numerical  simulation  of  the  directional 
solidification  process  with  syceep  temperature  gradient  were 
investigated.  In  order  to  improve  the  accuracy  of  calculation, 
the  effect  of  the  term  in  the  heat  transfer 

equation  on  the  outcome  with  steep  temperature  gradient  is 
discussed.  Problems  such  as  how  to  treat  latent  heat  of 
solidification,  thermal  properties  at  the  casting/mold  interface 
and  the  initial  condition  are  also  analyzed.  Both  finite 
difference  and  finite  element  methods  were  used  to  calculate  the 
same  casting  under  identical  initial  and  boundary  conditions  with 


the  same^ physical  parameters. 

Table  of  Symbols 


The  results  were  compared.* 

A 


T  temperature  (°C) 

t  time  (s) 

c  specific  heat  (cal/g°C) 

3 

p  density  (g/cm  ) 

k  conductivity  coefficient  (cal/s.cm°C) 

w  heat  released  by  internal  source  per  unit  volume  and  time 

3 

(cal/s.cm  ) 

j  latent  heat  (cal/g) 

fg  ratio  of  solidified  phase  to  total  amount 

Cj  =  -J  9fg/aT 

T^  liquid  phase  curve  temperature  (°C) 

T  solid  phase  curve  temperature  (®C) 

s  2  1/2 
b  thermal  storage  coefficient  (cal/cm  .s  .8C) 

Cp  physical  specific  heat  (cal/g. °C) 
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M  =  Ax.k/c.p. (Ax)  ,  At  is  the  time  interval  and  Ax  is  the 

step  length  in  space 

p  superscript  of  T  representing  time  interval 

i  step  number  in  z-direction 

j  step  number  in  r-direction 

m  slope  of  liquid  phase  curve 


Cq  percentage  of  solute  concentration  '  /31-4 

Icq  distribution  index 

Subscripts 

m  or  M  castings 
s  mold 

if  interface 

I.  Introduction 

Computer  numerical  simulation  of  the  solidification  process 
is  a  new  subject  of  study  in  casting.  In  recent  years,  its  scope 
and  width  have  been  expanding  rapidly.  One  of  the  directions  is 
to  study  the  directional  solidification  process  of  castings  by 
numerical  simulation. 

One  of  the  key  problems  in  conducting  numerical  simulation 
for  the  solification  process  is  the  selection  of  thermodynamic 
parameters  for  the  castings  and  molds.  After  taking  into  account 
that  there  is  very  little  data  on  directionally  solidified  high 
temperature  alloys  and  molds,  A1-4.5Cu  alloy  was  chosen  in  this 
work.  The  mold  is  also  metallic,  which  is  a  new  trend  in 
directional  solidification. 

In  the  solidification  process,  the  bottom  of  the  mold  was 
cooled  while  it  was  heated  all  the  way  around.  The  simulation 
condition  is  casting  temperature  740°C,  mold  preheating 
temperature  (measured  4cm  from  the  bottom)  650°C  and  water  cooled 
crystallized  (200  1/hr).  The  axial  temperature  gradient  was 
found  to  be  30-60°C/cm  experimentally.  The  specimen  is  a  3cm 
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diameter  10cm  long  cylinder.  The  columnar  structure  is  shown  in 
Figure  1[1]. 


Figure  1.  Columnar  Structure  of  the  Test  Ingot 

II.  Basic  Equations 

Due  to  the  fact  that  the  melting  point  of  the  aluminum  alloy 
used  in  the  experiment  is  not  very  high  (not  exceeding  750°C), 
the  temperature  differential  in  directional  solidification  is 
relatively  small.  After  neglecting  radiation  and  converting 
convection  to  conduction,  the  complicated  heat  transfer  problem 
is  simplified  to  a  simple  conduction  problem. 

The  basic  heat  transfer  equation  is  as  follows: 

e P-^--V(fcvr) 


(1) 


For  an  axial  symmetric  object,  equation  (1)  can  be  expressed 
in  cylindrical  coordinates 


»T  _  i  a  f  t  aT  \  ,  a  (  .  aT  \ 

Tr  r  -5rir*-s-J+-n-(‘-»r) 


Both  finite  difference  and  finite  element  methods  were  used 
to  calculate  the  same  casting  under  the  same  condition.  The 
aieshes  for  the  casting  and  the  mold  are  shown  in  Figures  2  and  3. 
The  discrete  elements  are  not  the  same  for  those  two  methods. 

From  equation (2) , the  numerical  equation  for  the  finite 
difference  method  is 


rft-  4MTt.l.l  +  MTt.ul+MTt.u,+  (l  “  (3) 

j  >  1 

T » "  ( 1 ' 20  -  1  +  [  1  +TU1-l))MTtn+MT?.u, 

+MTU,+  (l  -  4  M)Tl, 

0**7p"  7T*  (4) 


where  ek  *  at/cp  ak/aT. 
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Figure  2.  Mesh  for  Finite  Difference  Method 

1 .  casting 
__  2.  interface 

3.  mold 


Figure  3.  Mesh  for  Finite  Element  Method 

1  •  casting  2.  interface  3 


The  term  j-T£+.j  j/2ar)  in  equations  (3)  and  (4)  is  a  / 3 1 6 

result  of  taking  9k/9T(9T/9Z)2  into  consideration.  Because  the 
rate  of  change  of  k  with  respect  to  T  is  usually  small,  and  the 
temperature  gradient  in  sand  mold  casting  is  not  very  high,  the 
square  term  of  this  temperature  gradient  can  be  omitted  in 
numerical  calculation  using  finite  difference.  This  has  very 
little  effect  on  the  accuracy  of  calculation  for  sand  mold 
casting. 

However,  it  is  not  appropriate  to  neglect  the  temperature 
gradient  square  term  in  the  basic  equation  in  a  relatively  slow 
solidification  process  with  large  temperature  gradient.  The 
axial  temperature  gradient  is  always  very  large  in  casting 
ingots.  Whether  to  include  the  term  9k/9T(9T/9Z)2  has  a 
significant  effect  on  the  accuracy  of  the  calculation.  According 
to  the  results  shown  in  Figure  4,  the  calculated  result  deviates 
positively  from  the  experimental  value  when  0^  =  0,  i.e., 
neglecting  the  square  of  the  temperature  gradient  term,^'k  =  ~°o00002 
chosen  properly,  a  negative  deviation  would  result,  such  as  in 
the  case  =  -0.000094.  It  was  estimated  that  0^  *  -0.00002 
from  the  k  vs  T  curve  for  this  aluminum  alloy.  Strictly 
speaking,  should  also  be  a  function  of  T  because  k  is  a 
function  of  T. 

The  numerical  equation  for  equation  (2)  in  the  finite 
element  method: 

+  +  ^p]{n,.4,+  <QK.A,+  {Q>,-0cO 

where  H  is  a  matrix  related  to  the  thermal  conductivity 
coefficient,  P  is  a  matrix  related  to  specific  heat  and  density 
and  Q  is  series  related  to  the  natural  boundary  condition. 

Hence,  I T > T  can  be  obtained  by  solving  the  set  of  linear 
equations  when  {T}  .  is  known. 

t  —  a  T 

In  the  finite  element  method,  it  was  assumed  that  there  was 
no  heat  flux  across  axisymmetric  lines,  i.e.,  an  adiabatic 
process.  The  treatment  involves  the  use  of  a  third  type  of  boundary 
condition:  -k  9T/9r|axial  =  0^^.  In  our  work,  the  first  type 


of  boundary  condition  is  used  on  three  other  sides.  Along  the 
axis,  v  =  0.  Thus,  equation  (5)  can  be  simplified: 
r<  Df 
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Figure  4.  Comparison  Between  Simulated  Temperatures  with 

Different  Values  of  0V  and  its  Experimental  Data 
—  (at  Z  =  5.1cm) 

key,  1.  experimental  value 

[h  +  -^-p]iT).+  [h-~-t  p]<T>,.^=  0  (6)  / 3 1 7 

This  is  the  equation  used  in  this  work. 

Because  physical  parameters  of  the  material  such  as  k,  c, 

and  p  are  functions  of  T,  therefore,  the  matrices  [H]  and  [P] 

T41 

vary  with  T  . 

The  following  is  a  comparison  of  these  two  methods. 

The  explicit  finite  difference  method  is  more  accurate  when 
its  stability  condition  is  met.  Furthermore,  there  is  no  matrix 
in  the  process.  The  storage  requirement  is  small.  However,  the 
time  interval  is  limited  by  the  stability  requirement. 

Because  an  element  can  be  divided  arbitrarily  in  the  finite 
element  method,  it  is  more  superior  in  dealing  with  a  castings  of 
complex  shape.  A  simple  triangular  element  method  was  used  with 
linear  intrapolation,  therefore,  the  error  is  larger. 


Figure  5  is  a  comparison  of  these  two  methods  against  the 
experimental  result. 


Figure  5.  Comparison  Between  the  Results  of  the  Two 
Computation  Methods  and  Experimental  Data 

1  •  finite  difference  method 

2.  finite  element  method 

3.  experimental  value 

III.  Treatment  of  Several  Specific  Problems 


1 .  Latent  Heat 

Latent  heat  must  be  considered  in  the  heat  transfer  process 
of  cast  solidification.  The  latent  heat  of  crystallization  is 
high  for  metals  and  alloys.  For  example,  it  is  92.0  cal/g  for 
Al-4.5Cu.  If  the  liquid-solid  line  spacing  is  103°C  and  the 

31 . 9cal.  /  S  , 

specific  heat  is  0.31  cal/g°C,  then  the  heat  of  solidification^ 
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stage  will  affect  the  accuracy  of  calculation. 

If  there  is  a  heat  source  in  the  object,  then  the  expression 
for  equation  (1)  becomes 

cP-^«  V-(*VT)+u> 

*T  (7) 

Let  us  assume  that  the  latent  heat  .is  a  heat  source  in  the 
solidification  stage  and  the  amount  of  heat  released  is 
proportional  to  the  percentage  of  solid  precipitated' per  unit 
time,  i.e. , 

*T  (8) 

Usually,  f  is  also  a  function  of  T  and  T  is  a  function  of 
coordinate  and  time.  Therefore: 

ifs  _  ifs  iT 

dX  ~  »T  it 

—o'# 


Substituting  the  above  into  equation  (7),  we  get 


(c-/Sf)p4r-v<*vr> 


Let  c,~-  J 

(10) 

then,  the  specific  heat  c  in  the  solidification  process  should 
include  the  physical  specific  heat  Cp  and  the  latent  heat 
specific  heat  c j ,  i.e., 


c  «c,+c. 


(11) 


This  is  the  so-called  reduced  specific  heat  method  which  is 
applicable  to  alloys  of  a  certain  crystallization  temperature 


range 


The  simplest  treatment  is 


The  result  is  Cj  =  J/T^-Tg  which  is  a  mean  formula. 

If  we  combine  this  with  the  binary  diagram  of  an  alloy, 
another  type  of  formula  can  be  obtained.  For  instance, 


Tl-T 


cu 


«c.+rt-n 


(13) 


where  the  lever  law  in  equilibrium  was  used  in  the  derivation  of 
aquation  (13). 

The  distribution  of  the  solidified  solid  phase  is  related  to 
the  heat  exchange  of  the  system.  When  the  temperature  drops 
slowly  near  equilibrium,  existing  diagrams  can  be  used  to  obtain 
an  expression  for  fg.  The  situation  is  more  complicated  for  non¬ 
equilibrium  conditions. 

The  non-equilibrium  solidification  problem  in  directional 
crystallization  can  be  treated  according  to  the  normal 
solidification  equation 


or 


CL-CJlr‘ 
(Tm-T  V-L 


.to 


(14) 


where  Tm  is  the  melting  point  of  the  base  metal  in  the  alloy. 
The  relation  between  ^  and  temperature  T  can  be  obtained  based 
on 

Three  types  of  reduced  specific  heat  formulas  were  used  to 
simulate  the  same  casting  condition. 

(1)  reciprocal  square  formula  obtained  from  eq.  (14) 

mCJ _ 

Cj~  ( 1  — fe,)(mC,+7’t— T )* 

(2)  equation  of  average  derived  from  eq.  (13) 

Cj-  Tl-T, 

(3)  trapezoidal  equation,  i.e.,  the  slope  of  Cj  vs  T  line 
is  different  in  various  temperature  intervals. 


Equation  (10)  shows  that  the  simulation  of  latent  heat  is 
only  accurate  when  the  fg  vs  T  relation  is  known  for  the  entire 
solidification  process.  A  feasible  method  is  to  find  the  fg  vs  T 
curve  experimentally.  The  trapezoidal  equation  is  obtained  from 
the  fg-T  curve  of  the  Al-4.5Cu  alloy. 

Normally,  although  the  temperature  field  may  differ  by  using 
a  different  formula  to  simulate  the  latent  heat  of  the  case, 
however,  it  has  little  effect  on  the  outcome.  In  directional 
solidification,  external  conditions  are  controlled  rigorously. 

One  can  see  from  the  simulation  results  that  different  latent 
heat  formulas  would  yield  different  isothermal  solid-liquid 
interface  curvature  and  solidification  zone  width.  Therefore, 
the  solidification  of  latent  heat  must  be  handled  well  in 
directional  solidification.  Otherwise,  wrong  results  may  be 
obtained. 

Figure  6  shows  the  isothermal  surfaces  between  liquidus  and 
solidus  of  ingot  40  and  60  seconds  after  casting  with  different 
latent  haat  formulas.  The  *  symbol  respresents  an  actual 
experimental  point  of  liquidus  in  the  casting  and  the  symbol  • 
represents  an  experimental  solidus  point.  One  can  see  that  it  is 
more  appropriate  to  use  the  trapezoidal  formula. 

2.  Heat  Transfer  at  the  Interface  / 3 1 9 

For  convenience,  it  is  usually  classified  as  ideal  and  non¬ 
ideal  contacts.  In  reality,  ideal  contact  is  rare.  In  order  to 
improve  the  accuracy  of  calculation,  the  study  of  non-ideal 
contact  is  receiving  more  attention.  At  high  temperature,  the 
physico-chemical  process  at  the  interface  is  very  complex.  The 
mechanism  causing  heat  transfer  variation  is  not  very  clear. 

Mold  coating,  oxide  layer  on  the  case  and  gap  created  in 
solidification  will  change  the  heat  transfer  at  the  interface. 
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t«40iec 


4.  t>«0*k 


6. 

*«&X 


Figure  i>. 

Influence  of  Different  Formulae  for  Latent  Heat  on 

the 

Isothermal  Surface  Between  Liquidus  and  Solidus 

lUp 

1. 

liquidus 

2. 

solidified  phase 

3. 

solidus 

4. 

reciprocal  square  formula 

5. 

average  formula 

6. 

trapezoidal  formula 

Interfacial  heat  transfer  may  be  simulated  using  different 
formulas.  In  the  finite  difference  method,  a  reduced  thermal 
conductivity  coefficient  method  is  usually  used.  Its  formula  is 
as  follows: 


where  k  is  a  correction  factor. 


fci- 


rj-  2  kmkt 


(15) 

In  ideal  contact  situations,  k 


In  the  finite  element  method,  the  non-ideal  contact  is 
corrected  by  modifying  the  thermal  conductivity  coefficient  of 
the  mold  element  at  the  interface^ 

3.  Initial  Condition 

The  difficulty  in  determining  the  initial  condition  is  to 
calculate  the  temperature  at  the  mold  and  casting  interface.  In 
sand  casting,  the  following  formula  is  used  to  obtain  the  initial 
interfacial  temperature. 

t  _  bmTm,+btT„ 

T,t — Vu+br 

(16) 

where  tmi  and  are  the  initial  temperature  of  the  casting  and 
mold,  respectively. 

There  is  some  deviation  when  equation  (16)  is  used  in  a 
directional  solidification  process  with  steep  temperature 
gradient.  The  reason  is  that  the  formula  does  not  include  the 
latent  heat  released  by  some  solidified  metal  near  the  initial 
interface.  Equation  (16)  is  corrected  by  this  effect: 


where  b^  =  ^kMPM^cM  +  J^TL-Ts^ 
released  latent  heat  include 


«r  _  b’JT m, + btT „ 

"  ii+6, 

l 

.  Other  formulas  taking 


(17) 


t.—  P"c"7'mi  +  PsCsTtt+TLpuJ/(TL—Ts 
*  if  ^  _  n — F77m - m  v  1 


P*c*+P*c*+Pi»//(7't— Ts) 


(18) 


The  following  table  lists  the  initial  temperature  values 
calculated  with  three  above  formulas,  as  well  as  the  experimental 
value  (condition:  casting  temperature  of  the  ingot  and  pre¬ 
heating  temperature  of  the  crystallizer  are  known).  Equation 
(17)  gives  the  most  close  result  as  compared  to  the  measured 
value. 
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Temperature  °C 
formula 


condition 


casting  temperature  740 °C 
crystallizer  temperature  520°C 


Calculated  based  on 
Calculated  based  on 
Calculated  based  on 
Experimental 


eq 

(16) 

598 

eq 

(17) 

635 

eq 

(18) 

642 

630 

Conclusions 

1.  In  order  to  improve  the  accuracy  of  numerical  simulation  of 
the  solidification  process  with  steep  temperature  gradient,  the 
temperature  graient  square  terms  such  as  ak/8T(aT/aZ)*  in  the 
heat  transfer  equation  cannot  be  neglected. 

2.  Different  formulas  to  deal  with  latent  heat  will  have 
different  effects  on  the  isothermal  surface  and  zone  width. 

3.  The'calculation  of  the  initial  temperature  at  the  casting 
mold  interface  must  take  the  latent  heat  released  by  the  casting 
into  consideration  when  the  temperature  gradient  is  steep. 

4.  There  are  advantages  and  disadvantages  in  both  methods. 

Finite  difference  is  more  accurate  when  the  stability  requirement 
is  met.  It  is  suited  for  small  computers  because  the  storage 
requirement  is  small.  Finite  element  is  suited  for  castings  of 
complicated  shape  because  elements  can  be  divided  arbitrarily. 

References 

[1]  He  Qingfeng,  Effect  of  Solidification  Condition  on 
Structure  and  Properties  of  Al-4.5Cu  Alloy,  Master's 
Thesis  of  Northwestern  Polytechnical  University,  (1981). 

C2]  Zhu  Bofang,  Principle  and  Application  of  Finite  Element 
Method,  Hydroelectric  Power  Publishing  Co.,  (1979), 
pp  251-260. 

[33  Mathew,  J.  and  Brody,  H.D.,  Analysis  of  Heat  Transfer 
in  Continuous  Casting  Using  a  Finite  Element  Method, 


I 

Conference  on  Computer  Simulation  for  Materials  Appli¬ 
cations,  (1976),  p.  1142. 

;  C4]  Comini,  G.,  Delguidice,  S.,  Lewis,  R.W.  and  Zienkiewicz, 

;  O.Z.,  Finite  Element  Solution  of  Non-Linear  Heat  Conduction 

I  Problems  with  Special  References  to  Phase  Change,  Inter. 

|  Journal  of  Numerical  Methods  in  Eng.,  Vol.  8,  (1974), 

pp.  613-624. 

[5]  Guo  Keren,  Jin  Junzhe,  et  al,  Mathematical  Modeling  of 
Solidification  of  Large  Scale  Casting,  Dalian  Engineering 

|  Institute  Journal,  2(1980),  pp  1-17. 

[6]  Clyne,  T.W.  and  Davies,  G.J.,  The  Influence  of  Composition 
on  Solidification  Cracking  Susceptibility  in  Binary  Alloy 
Systems,  The  British  Foundreyman,  Vol.  74,  Apr.  (1981), 

l  p.  65. 

[7]  Pehlke,  R.D.,  Kirt,  M.J.,  Marrone,  R.E.,  and  Cook,  D.J., 
Numerical  Simulation  of  Casting  Solidification,  AFS  Cast 
Metals  Research  Journal,  Vol.  9,  June  (1973),  pp  49-55. 


15 


